%DASOV24.m

% Data analysis 

clear
tic

% GOLD/SILVER/PLATINUM/COPPER/NICKEL  - monthly averages from daily data
% Source: Bloomberg Daily data from London, 
%               except for Platinum use XPT for longer sample

% create cell-array Xall monthly (Xalld  daily)
% {1}   {2}     {3}     {4}     {5}
%  Gold Silver Platinum Copper Nickel
%  Each cell contains 3 times T matrix
%           1     2     3      
%       [ date  price price_inflation_adjusted ]


X =          xlsread('GoldBB.xlsx', 1, 'B9:B13455');
[~, timlab]= xlsread('GoldBB.xlsx',1,'A9:A13455');
Xalld{1} = flip([ datenum(timlab,'mm/dd/yyyy') X]); 

longsamgold=1;     % longsamgold=0 -- common sample length for all metals 87.2 (selected below)
                   % longsamgold=1 -- longest sample only gold 75/79

X =          xlsread('SilverBB.xlsx',1,'B9:B13645');
[~, timlab]= xlsread('SilverBB.xlsx',1,'A9:A13645');
Xalld{2} = flip([ datenum(timlab,'mm/dd/yyyy') X]); 

X =          xlsread('PlatinumBBXPT.xlsx',1,'B61:B9748');
[~, timlab]= xlsread('PlatinumBBXPT.xlsx',1,'A61:A9748');
Xalld{3} = flip([ datenum(timlab,'mm/dd/yyyy') X]); 

X =         xlsread('CopperBB.xlsx',1,'B9:B9573');
[~, timlab]= xlsread('CopperBB.xlsx',1,'A9:A9573');
Xalld{4} = flip([ datenum(timlab,'mm/dd/yyyy') X]); 

X =         xlsread('NickelBB.xlsx',1,'B9:B9376');
[~, timlab]= xlsread('NickelBB.xlsx',1,'A9:A9376');
Xalld{5} =flip([ datenum(timlab,'mm/dd/yyyy') X]); 

% Convert to monthly
for k=1:5
    Xm =  ConvDaytoMonth(Xalld{k});
    Xall{k}=Xm;
end

% Consumer price index
CPIr            = xlsread('CPIAUCSL.xls', 1, 'B12:B936');
[~, timlab]    = xlsread('CPIAUCSL.xls', 1, 'A12:A936');
dCPIr = datenum(timlab,'mm/dd/yyyy');

save cpidat CPIr dCPIr

% Deflate metal prices by CPI
for k=1:5
    N = max(size(Xall{k}))-1;     
    Xall{k} = [Xall{k} ones(N+1,1) ];
    for j=1:N
        Xall{k}(j,3)= Xall{k}(j,2)/CPIr( dCPIr ==Xall{k}(j,1));
    end
    Xall{k}(end,:)=[];    
end

% Load yields and inflation expectations

% Treasury Rates: 1 and 10 years
Y1r            = xlsread('GS1.xls', 1, 'B12:B862');
[~, timlab]    = xlsread('GS1.xls', 1, 'A12:A862');
dY1r = datenum(timlab,'mm/dd/yyyy');    
    
Y10r           = xlsread('GS10.xls', 1, 'B12:B862');
[~, timlab]    = xlsread('GS10.xls', 1, 'A12:A862');
dY10r = datenum(timlab,'mm/dd/yyyy'); 

Y10TIPSr       = xlsread('FII10.xls', 1, 'B12:B265');
[~, timlab]    = xlsread('FII10.xls', 1, 'A12:A265');
dY10TIPSr = datenum(timlab,'mm/dd/yyyy'); 
    
% One year inflation forecast from SPF: 75-80 GDPdef, after CPI
% Quarterly data 75Q1 to 24Q1
CPI1er           = xlsread('inflation24.xlsx', 1, 'f22:f218');
Nm=    12*(2024-1975)+3;

CPI1e = zeros(Nm,1);
% Assign to in quarter months 
CPI1e(1:3:Nm-2)  = CPI1er;
CPI1e(2:3:Nm-1)  = CPI1er;
CPI1e(3:3:Nm)    = CPI1er;

dCPI1e = datenum(1975, 1 + [1:Nm]' - 1, 1); % First day of the month


% 10-Year Inflation expectation starting SPF 1991.Q4 -- 1991.10 unit 2024.q1
CPI10er           = xlsread('inflation24.xlsx', 1, 'e89:e218');
Nm = 12*(2024-1991)-6;
CPI10e = zeros(Nm,1);
% Assign to in quarter months
CPI10e(1:3:Nm-2)    = CPI10er;
CPI10e(2:3:Nm-1)    = CPI10er;
CPI10e(3:3:Nm)      = CPI10er;

% Inflation expectations extensions 1979.4 - 1991.3 Blue Chip: March and October
% 2 versions to deal with having only 2 annual data points:
% (1)Column'e' keep values until next available, (2) 'f' average (future info)
CPI10Xer           = xlsread('additional-CPIE10.xlsx', 1, 'e15:e62');
Nmex = 12*(1991-1979);
CPI10eX=zeros(Nmex,1);
% Assign to in quarter months
CPI10eX(1:3:Nmex-2)    = CPI10Xer;
CPI10eX(2:3:Nmex-1)    = CPI10Xer;
CPI10eX(3:3:Nmex)      = CPI10Xer;

CPI10e = [CPI10eX; CPI10e];
dCPI10e = datenum(1979, 10 + [1:(Nm+Nmex)]' - 1, 1); % First day of the month

% Inflation adjust yields
% 1Y from 1975 Jan to 2023 Dec
dY1IA = datenum(1975, 1 + [1:(12*(2023-1975)+12)]' - 1, 1); % First day of the month

Y1IA  = 100*((1+Y1r(dY1r>=dY1IA(1)&dY1r<=dY1IA(end))/100)./...
                    (1+ CPI1e(dCPI1e>=dY1IA(1)&dCPI1e<=dY1IA(end)) /100)-1) ;
% 
% 10Y from 1979.10 to 2023.12
dY10IA = datenum(1979, 10 + [1:(12*(2023-1979)+3)]' - 1, 1); % First day of the month
Y10IA  = 100*((1+Y10r(dY10r>=dY10IA(1)&dY10r<=dY10IA(end))/100)./...
                    (1+ CPI10e(dCPI10e>=dY10IA(1)&dCPI10e<=dY10IA(end)) /100)-1) ;

% Replace by TIPS: 2003.1 to 2023.12
Y10IA(dY10IA>=datenum(2003,1,1) & dY10IA<=datenum(2023,12,1)) = ...
    Y10TIPSr(dY10TIPSr>=datenum(2003,1,1) &dY10TIPSr<=datenum(2023,12,1)); 


figure(1)
p=plot(Xall{1}(:,1), 2.5*Xall{1}(:,3),...
       dY10IA, Y10IA,...
       dY1IA, Y1IA );
dateFormat = 12;
datetick('x',dateFormat)
p(1).Color = [184, 134, 11]/255; %Dark goldenrod / #b8860b hex color
p(1).LineWidth = 1.5;
p(2).Color = [13, 71, 161]/255; % #0D47A1
p(2).LineWidth = 1.5;
p(3).Color = [79, 195, 247]/255;
p(3).LineWidth = 1.5;
title(['Real Gold Price and Real Treasury Rates & TIPS ' ])
axis([datenum(1975,1,1) datenum(2023,12,31) -5 25])
legend('Gold', '10-year Treasury/TIPS','1-year Treasury','Location','NorthEast') 
grid

% exportgraphics(gcf,'fig1.png','Resolution',300)

% Regressions of log price change on log yield change conditional on low/high yields
% 1Y and 10Y here, and inside Runregsf(.):
% 1) Annual and monthly
% 2) 33% and 50% cutoff

% Select a common sample

sambeg = datenum(1987,2,1);    %XPT
samend = datenum(2023,12,1);

% *** Gold ***
% *** 1Y rate ****
pric = Xall{1}(:,3);
dpric = Xall{1}(1:end,1);
yiel = Y1IA(1:end);
dyiel = dY1IA(1:end);
TGold1Y = Runregsf(pric,yiel,dpric,dyiel,sambeg,samend);
% *** 10Y rate ****
yiel = Y10IA(1:end);
dyiel = dY10IA(1:end);
TGold10Y = Runregsf(pric,yiel,dpric,dyiel,sambeg,samend);

if longsamgold==1
    % *** Gold ***
    % *** 1Y rate ****
    pric = Xall{1}(:,3);
    dpric = Xall{1}(1:end,1);
    yiel = Y1IA(1:end);
    dyiel = dY1IA(1:end);
    % Select sample manually -- largest possible for monthly change
    sambeggo = datenum(1975,1,1);
    samendgo = datenum(2023,12,1);
    TGold1Y = Runregsf(pric,yiel,dpric,dyiel,sambeggo,samendgo);
    % *** 10Y rate ****
    yiel = Y10IA(1:end);
    dyiel = dY10IA(1:end);
    % Select sample manually -- largest possible for monthly change
    sambeggo = datenum(1979,10,1);
    samendgo = datenum(2023,12,1);
    TGold10Y = Runregsf(pric,yiel,dpric,dyiel,sambeggo,samendgo);
end


% *** Silver ***
% *** 1Y rate ****
pric = Xall{2}(:,3);
dpric = Xall{2}(1:end,1);
yiel = Y1IA(1:end);
dyiel = dY1IA(1:end);
TSilver1Y = Runregsf(pric,yiel,dpric,dyiel,sambeg,samend);
% *** 10Y rate ****
yiel = Y10IA(1:end);
dyiel = dY10IA(1:end);
TSilver10Y = Runregsf(pric,yiel,dpric,dyiel,sambeg,samend);

% *** Platinum ***
% *** 1Y rate ****
pric = Xall{3}(:,3);
dpric = Xall{3}(1:end,1);
yiel = Y1IA(1:end);
dyiel = dY1IA(1:end);
% Select sample manually -- largest possible for monthly change
TPlatinum1Y = Runregsf(pric,yiel,dpric,dyiel,sambeg,samend);
% *** 1Y rate ****
yiel = Y10IA(1:end);
dyiel = dY10IA(1:end);
TPlatinum10Y = Runregsf(pric,yiel,dpric,dyiel,sambeg,samend);

% *** Copper ***
% *** 1Y rate ****
pric = Xall{4}(:,3);
dpric = Xall{4}(1:end,1);
yiel = Y1IA(1:end);
dyiel = dY1IA(1:end);
TCopper1Y = Runregsf(pric,yiel,dpric,dyiel,sambeg,samend);
% *** 1Y rate ****
yiel = Y10IA(1:end);
dyiel = dY10IA(1:end);
TCopper10Y = Runregsf(pric,yiel,dpric,dyiel,sambeg,samend);

% *** Nickel ***
% *** 1Y rate ****
pric = Xall{5}(:,3);
dpric = Xall{5}(1:end,1);
yiel = Y1IA(1:end);
dyiel = dY1IA(1:end);
TNickel1Y = Runregsf(pric,yiel,dpric,dyiel,sambeg,samend);
% *** 1Y rate ****
yiel = Y10IA(1:end);
dyiel = dY10IA(1:end);
TNickel10Y = Runregsf(pric,yiel,dpric,dyiel,sambeg,samend);


TGold1Y
TGold10Y
if longsamgold==1
    'for longer gold sample'
end

TSilver1Y
TSilver10Y

TPlatinum1Y
TPlatinum10Y

TCopper1Y
TCopper10Y

TNickel1Y
TNickel10Y


% Produce plots of rsquares and slopes for different medals
% plotsfivemetals

plotsfivemetals50
% exportgraphics(gcf,'fig4.png','Resolution',300)

% Gold prices function of yields scatterplots with fitted lines
% 1 year rates
sambeggo1Y = datenum(1975,1,1);
samendgo1Y = datenum(2023,12,1);
GPIA = Xall{1}(Xall{1}(:,1)>=sambeggo1Y&Xall{1}(:,1)<=samendgo1Y,3) ;

% 1 is 75.1
sampf = ones(numel(GPIA),1);
cuto1 = 0;      
samp1 = (Y1IA<cuto1);

stats = regstats2(GPIA(sampf==1),[Y1IA(sampf==1)],'linear');                                    yh1 = stats.yhat;
stats = regstats2(GPIA(sampf==1),[Y1IA(sampf==1) (1-samp1).*(Y1IA(sampf==1)-cuto1)],'linear');  yh2 = stats.yhat;  rsq2=stats.adjrsquare;
stats = regstats2(GPIA(sampf==1),[Y1IA(sampf==1) Y1IA(sampf==1).^2 Y1IA(sampf==1).^3 Y1IA(sampf==1).^4  Y1IA(sampf==1).^5 Y1IA(sampf==1).^6],'linear');   beta3=stats.beta; pval3=stats.hac.pval; rsq3=stats.adjrsquare; yh3 = stats.yhat;
stats = regstats2(GPIA(sampf==1),[Y1IA(sampf==1) Y1IA(sampf==1).^2 Y1IA(sampf==1).^3 Y1IA(sampf==1).^4  ],'linear');   beta3=stats.beta; pval3=stats.hac.pval; rsq3=stats.adjrsquare; yh3 = stats.yhat;

figure(21)
plot(Y1IA(sampf==1), GPIA(sampf==1),'.',...
Y1IA(sampf==1),yh1 ,'-',...
Y1IA(sampf==1),yh2 ,'.',...
Y1IA(sampf==1),yh3 ,'.' )
title('Real Gold Price as a Function of 1-year Real Rate')

% exportgraphics(gcf,'fig2.png','Resolution',300)


% 10 year rates
sambeggo10Y = datenum(1979,10,1);
samendgo10Y = datenum(2023,12,1);
GPIA10 = Xall{1}(Xall{1}(:,1)>=sambeggo10Y&Xall{1}(:,1)<=samendgo10Y,3) ;

samp10f = ones(numel(GPIA10),1);
cuto10 = 2.2 ;     
samp10 = (Y10IA<cuto10);

stats = regstats2(GPIA10(samp10f==1),[Y10IA(samp10f==1)],'linear');  yh1 = stats.yhat;
stats = regstats2(GPIA10(samp10f==1),[Y10IA(samp10f==1) (1-samp10).*(Y10IA(samp10f==1)-cuto10)],'linear');   yh2 = stats.yhat; rsq2=stats.adjrsquare;
stats =  regstats2(GPIA10(samp10f==1),[Y10IA(samp10f==1) Y10IA(samp10f==1).^2 Y10IA(samp10f==1).^3 Y10IA(samp10f==1).^4  Y10IA(samp10f==1).^5 Y10IA(samp10f==1).^6],'linear'); beta3=stats.beta; pval3=stats.hac.pval; rsq3=stats.adjrsquare; yh3 = stats.yhat;
stats =  regstats2(GPIA10(samp10f==1),[Y10IA(samp10f==1) Y10IA(samp10f==1).^2 Y10IA(samp10f==1).^3 Y10IA(samp10f==1).^4  ],'linear'); beta3=stats.beta; pval3=stats.hac.pval; rsq3=stats.adjrsquare; yh3 = stats.yhat;

figure(22)
plot(Y10IA(samp10f==1), GPIA10(samp10f==1),'.',...
Y10IA(samp10f==1),yh1 ,'-',...
Y10IA(samp10f==1),yh2 ,'.',...
Y10IA(samp10f==1),yh3 ,'.' )
title('Real Gold Price as a Function of 10-year Real Rate')

% exportgraphics(gcf,'fig3.png','Resolution',300)

% SAVED for model Simulations / Regressions
begm10 =dY10IA; 
Y1IA10=Y1IA(dY1IA>=dY10IA(1)&dY1IA<=dY10IA(end));  
save datforSuff24 GPIA10 Y10IA begm10 Y1IA10


Movday =        xlsread('MERmove.xlsx',1,'I10:K5365');
[~, timlab]=    xlsread('MERmove.xlsx',1,'H10:H5365');
Movday =     [ datenum(timlab,'mm/dd/yyyy') Movday]; 

% Convert to monthly
Movm =  ConvDaytoMonth(Movday);

% SAVED for model Simulations / Regressions
save movdat Movm


%% Futures-DAILY data will average Starting 4/02/1979 - 12/31/2023: 
Fuc1f2Dr             = xlsread('GoldFut.xlsx', 1, 'G306:G11571');
Fuc1f8Dr=Fuc1f2Dr;      % have two 1-month vectors comparing to 2 or 8 months
Fuc2Dr             = xlsread('GoldFut.xlsx', 1, 'H306:H11571');
Fuc8Dr             = xlsread('GoldFut.xlsx', 1, 'N306:N11571');
[~, timlab]      =   xlsread('GoldFut.xlsx', 1, 'F306:F11571');
dFu2Dr = datenum(timlab,'mm/dd/yyyy');
dFu8Dr = datenum(timlab,'mm/dd/yyyy');

% get rid of days with missing values
induf2 = zeros(max(size(Fuc1f2Dr)),1);
for j = 1: max(size(Fuc1f2Dr))
    if isnan(Fuc2Dr(j) -Fuc1f2Dr(j) )
        induf2(j)=1;
    end
end
Fuc1f2Dr(induf2==1)=[];
Fuc2Dr(induf2==1)=[];
dFu2Dr(induf2==1)=[];

induf8 = zeros(max(size(Fuc1f8Dr)),1);
for j = 1: max(size(Fuc1f8Dr))
    if isnan(Fuc8Dr(j) -Fuc1f8Dr(j) )
        induf8(j)=1;
    end
end
Fuc1f8Dr(induf8==1)=[];
Fuc8Dr(induf8==1)=[];
dFu8Dr(induf8==1)=[];

% vector first day of each month from 1979.4 to 2023.12: 44 years 9 month
begmfu = zeros(45*12,1);
j=1;
for yea = 1979:2023
    for mon = 1:12
        begmfu(j,1) =   datenum(yea,mon,1);
        j=j+1;
    end
end
begmfu(1:3)=[];
Tfu = 44*12+9;

begmfu(Tfu+1,1) =datenum(2024,1,1);

Fuc1f2 = zeros(Tfu,1);
Fuc1f8 = zeros(Tfu,1);
Fuc2 = zeros(Tfu,1);
Fuc8 = zeros(Tfu,1);

for j = 1:Tfu
    Fuc1f2(j,1) = mean( Fuc1f2Dr(dFu2Dr >= begmfu(j)&dFu2Dr < begmfu(j+1))  );
    Fuc1f8(j,1) = mean( Fuc1f8Dr(dFu8Dr >= begmfu(j)&dFu8Dr < begmfu(j+1))  );  
    Fuc2(j,1)  = mean( Fuc2Dr(dFu2Dr >= begmfu(j)&dFu2Dr < begmfu(j+1))  );
    Fuc8(j,1)  = mean( Fuc8Dr(dFu8Dr >= begmfu(j)&dFu8Dr < begmfu(j+1))  );
end

% Need 1 month T-bill Rate
TB4wksraw= xlsread('CRSP26wkTermStr.xlsx',1,'e1:e11244');
[~, timlab]      =   xlsread('CRSP26wkTermStr.xlsx',1,'b1:b11244');
dTB4wksraw = datenum(timlab,'mm/dd/yyyy');

TB4wks = zeros(Tfu,1);
for j = 1:Tfu
    TB4wks(j,1) = mean( TB4wksraw(dTB4wksraw >= begmfu(j)&dTB4wksraw < begmfu(j+1))  );
end

Furat1M10 = (Fuc2(7:Tfu)./Fuc1f2(7:Tfu))./((1+ TB4wks(7:Tfu)).^(1/12));

%Maturity adjustments for 1 year horizon
Futmatadj = ones(Tfu,1);        % Apr 8-1gap is 12 month, May 11month etc, Dec 12
Futmatadj(2:2:Tfu-1) = (12/11)*ones((Tfu-1)/2,1);
Furat1Y = (Futmatadj.*(Fuc8./Fuc1f8-1)+1)./ (1+ Y1r(dY1r >= datenum(1979,4,1)&dY1r <= datenum(2023,12,1))/100);
Furat1Y10 =   Furat1Y(7:Tfu);       % cut to 1979.10 (since 79.4)

% SAVED for model Simulations ; cut before to 1979.10 (since 79.4)
save datfutforSuff24MON Furat1M10 Furat1Y10


toc





